##  This file contains the code to replicate the analyses reported in Tables 3 and 4
##  and Figures 2 and 3 in the article
##

rm(list = ls(all.names = TRUE))
library(Hmisc)
library(lmtest)
library(sandwich)
library(pscl)
library(MASS)
library(mice)
library(miceadds)
library(openxlsx)


#setwd("~/Dropbox/Current projects/Nazi trials/analysis")
load("Nazi trials replication data.RData")
file.info("Nazi trials replication data.RData")$mtime


##  how many unique trials?
length(unique(cases.temp$ID))

##  how many unique defendants?
length(unique(cases.temp$unique.defendant.ID))

##  how many convictions?
table(cases.temp$convicted)
round(mean(cases.temp$convicted),2)


##  source helper functions
source("probit2.r")



##
##  Estimate models and save the results
##  Also conduct Monte Carlo simulations of first differences and plot the results
##  For now we save results in one big table; below this table gets split into 2 separate tables
tab1 <- matrix("",29,6)
colnames(tab1) <- paste0("(",1:ncol(tab1),")")

rownames(tab1) <- c(
					"% judges who were appointed in 1933--45", "",
					"% judges who were appointed pre-1933", "",
					"% judges who completed legal training in 1933--45", "",
					"% judges who completed legal training pre-1933", "",
					"avg. share of judges' legal education in 1933--45", "",
					"avg. share of judges' legal education pre-1933", "",
					"% prosecutors who were appointed in 1933--45", "",
					"% prosecutors who were appointed pre-1933", "",
					"% prosecutors who completed legal training in 1933--45", "",
					"% prosecutors who completed legal training pre-1933", "",
					"avg. share of prosecutors' legal education in 1933--45", "",
					"avg. share of prosecutors' legal education pre-1933", "",
					"Year FE","State FE", "Crime categorization FE",
					"$p$-value Wald test % judges","$p$-value Wald test % prosecutors")

sims.tab <- vector(mode = "list", length = 6)


##  add additional predictors?
add.pred <- " + state"	#  state FEs; all models include crime categories and year FEs



##  col 1
add.pred.MI <- c("NS_judge_@_20","NS_prosecutor_@_20")

##  probit model
out <- glm.MI(d = cases.temp, outcome = "convicted", add.pred.MI = add.pred.MI, add.pred = add.pred,
sims = 1, wald = NULL)
tab1[c(1,13),1] <- out[[1]][2:3,1]
tab1[c(2,14),1] <- out[[1]][2:3,2]
sims.tab[[1]][1] <- list(out[[3]])



##  col 2
add.pred.MI <- c("NS_judge_@_20","NS_prosecutor_@_20","NS_judge_@_19","NS_prosecutor_@_19")

##  probit
out <- glm.MI(d = cases.temp, outcome = "convicted", add.pred.MI = add.pred.MI, add.pred = add.pred,
sims = 1, wald = list(c(1,3),c(2,4)))
tab1[c(1,13,3,15),2] <- out[[1]][2:5,1]
tab1[c(2,14,4,16),2] <- out[[1]][2:5,2]
tab1[28,2] <- out[[2]][1]
tab1[29,2] <- out[[2]][2]
sims.tab[[2]][1] <- list(out[[3]])

out <- glm.MI(d = cases.temp, outcome = "convicted", add.pred.MI = add.pred.MI, add.pred = add.pred,
sims = 3, wald = list(c(1,3),c(2,4)))
sims.tab[[2]][2] <- list(out[[3]])



##  col 3
add.pred.MI <- c("NS_judge_@_2","NS_prosecutor_@_2")

##  probit
out <- glm.MI(d = cases.temp, outcome = "convicted", add.pred.MI = add.pred.MI, add.pred = add.pred,
sims = 1, wald = NULL)
tab1[c(5,17),3] <- out[[1]][2:3,1]
tab1[c(6,18),3] <- out[[1]][2:3,2]
sims.tab[[3]][1] <- list(out[[3]])



##  col 4
add.pred.MI <- c("NS_judge_@_2","NS_prosecutor_@_2","NS_judge_@_5","NS_prosecutor_@_5")

##  probit
out <- glm.MI(d = cases.temp, outcome = "convicted", add.pred.MI = add.pred.MI, add.pred = add.pred,
sims = 1, wald = list(c(1,3),c(2,4)))
tab1[c(5,17,7,19),4] <- out[[1]][2:5,1]
tab1[c(6,18,8,20),4] <- out[[1]][2:5,2]
tab1[28,4] <- out[[2]][1]
tab1[29,4] <- out[[2]][2]
sims.tab[[4]][1] <- list(out[[3]])

out <- glm.MI(d = cases.temp, outcome = "convicted", add.pred.MI = add.pred.MI, add.pred = add.pred,
sims = 3, wald = list(c(1,3),c(2,4)))
sims.tab[[4]][2] <- list(out[[3]])



##  col 5
add.pred.MI <- c("NS_judge_@_24","NS_prosecutor_@_24")

##  probit
out <- glm.MI(d = cases.temp, outcome = "convicted", add.pred.MI = add.pred.MI, add.pred = add.pred,
sims = 1, wald = NULL)
tab1[c(9,21),5] <- out[[1]][2:3,1]
tab1[c(10,22),5] <- out[[1]][2:3,2]
sims.tab[[5]][1] <- list(out[[3]])



##  col 6
add.pred.MI <- c("NS_judge_@_24","NS_prosecutor_@_24","NS_judge_@_23","NS_prosecutor_@_23")

##  probit
out <- glm.MI(d = cases.temp, outcome = "convicted", add.pred.MI = add.pred.MI, add.pred = add.pred,
sims = 1, wald = list(c(1,3),c(2,4)))
tab1[c(9,21,11,23),6] <- out[[1]][2:5,1]
tab1[c(10,22,12,24),6] <- out[[1]][2:5,2]
tab1[28,6] <- out[[2]][1]
tab1[29,6] <- out[[2]][2]
sims.tab[[6]][1] <- list(out[[3]])

out <- glm.MI(d = cases.temp, outcome = "convicted", add.pred.MI = add.pred.MI, add.pred = add.pred,
sims = 3, wald = list(c(1,3),c(2,4)))
sims.tab[[6]][2] <- list(out[[3]])



##  add checkmarks to table
tab1[25:27,1:6]	<- "\\checkmark"

##  split results table into two and omit irrelevant columns
tab1[c(1:4,13:16,25:29),1:2]					## TABLE 3 in article
tab1[c(9:12,21:24,25:29),5:6]					## TABLE 4 in article

latex(tab1[c(1:4,13:16,25:29),1:2], file = "")	## TABLE 3 in article
latex(tab1[c(9:12,21:24,25:29),5:6], file = "") ## TABLE 4 in article



##  plot first differences
sims.tab


#pdf(file = "plots/FD_appointed.pdf")
plot(c(1,1.925,2.075),c(
sims.tab[[1]][[1]][1],
sims.tab[[2]][[1]][1],
sims.tab[[2]][[2]][1]), xlab = "Model specification", ylab = "First difference in Pr(Conviction)",
main = "", ylim = c(-.25,.25), pch = c(19,19,19), bg = "black", cex.lab = 1.2,
xaxt = "n", xlim = c(0.75,2.325))

axis(side = 1, at = 1:2, labels = paste0("(",1:2,")"), cex.axis = 1.1)

abline(h = 0, col = "black", lty = 1, lwd = .7)

lines(c(1,1), c(sims.tab[[1]][[1]][1], sims.tab[[1]][[1]][1] + qnorm(.975)*sims.tab[[1]][[1]][2]))
lines(c(1,1), c(sims.tab[[1]][[1]][1], sims.tab[[1]][[1]][1] - qnorm(.975)*sims.tab[[1]][[1]][2]))


lines(c(1.925,1.925), c(sims.tab[[2]][[1]][1], sims.tab[[2]][[1]][1] + qnorm(.975)*sims.tab[[2]][[1]][2]))
lines(c(1.925,1.925), c(sims.tab[[2]][[1]][1], sims.tab[[2]][[1]][1] - qnorm(.975)*sims.tab[[2]][[1]][2]))

lines(c(2.075,2.075), c(sims.tab[[2]][[2]][1], sims.tab[[2]][[2]][1] + qnorm(.975)*sims.tab[[2]][[2]][2]))
lines(c(2.075,2.075), c(sims.tab[[2]][[2]][1], sims.tab[[2]][[2]][1] - qnorm(.975)*sims.tab[[2]][[2]][2]))

text(0.950, sims.tab[[1]][[1]][1],"% judges appointed in 1933-45", srt = 90, cex = .8)

text(1.875, sims.tab[[2]][[1]][1]-0.04,"% judges appointed in 1933-45", srt = 90, cex = .8)
text(2.125, sims.tab[[2]][[2]][1]+0.04,"% judges appointed pre-1933", srt = 90, cex = .8)

dev.off()



#pdf(file = "plots/FD_trained.pdf")
plot(c(1,1.925,2.075),c(
sims.tab[[5]][[1]][1],
sims.tab[[6]][[1]][1],
sims.tab[[6]][[2]][1]), xlab = "Model specification", ylab = "First difference in Pr(Conviction)",
main = "", ylim = c(-.35,.35), pch = c(19,19,19), bg = "black", cex.lab = 1.2,
xaxt = "n", xlim = c(0.75,2.325))

axis(side = 1, at = 1:2, labels = paste0("(",1:2,")"), cex.axis = 1.1)

abline(h = 0, col = "black", lty = 1, lwd = .7)


lines(c(1,1), c(sims.tab[[5]][[1]][1], sims.tab[[5]][[1]][1] + qnorm(.975)*sims.tab[[5]][[1]][2]))
lines(c(1,1), c(sims.tab[[5]][[1]][1], sims.tab[[5]][[1]][1] - qnorm(.975)*sims.tab[[5]][[1]][2]))


lines(c(1.925,1.925), c(sims.tab[[6]][[1]][1], sims.tab[[6]][[1]][1] + qnorm(.975)*sims.tab[[6]][[1]][2]))
lines(c(1.925,1.925), c(sims.tab[[6]][[1]][1], sims.tab[[6]][[1]][1] - qnorm(.975)*sims.tab[[6]][[1]][2]))

lines(c(2.075,2.075), c(sims.tab[[6]][[2]][1], sims.tab[[6]][[2]][1] + qnorm(.975)*sims.tab[[6]][[2]][2]))
lines(c(2.075,2.075), c(sims.tab[[6]][[2]][1], sims.tab[[6]][[2]][1] - qnorm(.975)*sims.tab[[6]][[2]][2]))

text(0.900, sims.tab[[5]][[1]][1],"avg. share of judges'",srt = 90, cex = .8)
text(0.950, sims.tab[[5]][[1]][1],"legal training in 1933-45",srt = 90, cex = .8)

text(1.825, sims.tab[[6]][[1]][1],"avg. share of judges'",srt = 90, cex = .8)
text(1.875, sims.tab[[6]][[1]][1],"legal training in 1933-45",srt = 90, cex = .8)

text(2.125, sims.tab[[6]][[2]][1]-.06,"avg. share of judges'",srt = 90, cex = .8)
text(2.175, sims.tab[[6]][[2]][1]-.06,"legal training pre-1933",srt = 90, cex = .8)

dev.off()



##  percentiles reported at bottom of Figure 3
# 1933-45 share
round(mean(
quantile(cases.temp$NS_judge_1_20,.25),
quantile(cases.temp$NS_judge_2_20,.25),
quantile(cases.temp$NS_judge_3_20,.25),
quantile(cases.temp$NS_judge_4_20,.25),
quantile(cases.temp$NS_judge_5_20,.25),
quantile(cases.temp$NS_judge_6_20,.25),
quantile(cases.temp$NS_judge_7_20,.25),
quantile(cases.temp$NS_judge_8_20,.25),
quantile(cases.temp$NS_judge_9_20,.25),
quantile(cases.temp$NS_judge_10_20,.25)),2)

round(mean(
quantile(cases.temp$NS_judge_1_20,.75),
quantile(cases.temp$NS_judge_2_20,.75),
quantile(cases.temp$NS_judge_3_20,.75),
quantile(cases.temp$NS_judge_4_20,.75),
quantile(cases.temp$NS_judge_5_20,.75),
quantile(cases.temp$NS_judge_6_20,.75),
quantile(cases.temp$NS_judge_7_20,.75),
quantile(cases.temp$NS_judge_8_20,.75),
quantile(cases.temp$NS_judge_9_20,.75),
quantile(cases.temp$NS_judge_10_20,.75)),2)


# pre-1933 share
# 1933-45 share
round(mean(
quantile(cases.temp$NS_judge_1_19,.25),
quantile(cases.temp$NS_judge_2_19,.25),
quantile(cases.temp$NS_judge_3_19,.25),
quantile(cases.temp$NS_judge_4_19,.25),
quantile(cases.temp$NS_judge_5_19,.25),
quantile(cases.temp$NS_judge_6_19,.25),
quantile(cases.temp$NS_judge_7_19,.25),
quantile(cases.temp$NS_judge_8_19,.25),
quantile(cases.temp$NS_judge_9_19,.25),
quantile(cases.temp$NS_judge_10_19,.25)),2)

round(mean(
quantile(cases.temp$NS_judge_1_19,.75),
quantile(cases.temp$NS_judge_2_19,.75),
quantile(cases.temp$NS_judge_3_19,.75),
quantile(cases.temp$NS_judge_4_19,.75),
quantile(cases.temp$NS_judge_5_19,.75),
quantile(cases.temp$NS_judge_6_19,.75),
quantile(cases.temp$NS_judge_7_19,.75),
quantile(cases.temp$NS_judge_8_19,.75),
quantile(cases.temp$NS_judge_9_19,.75),
quantile(cases.temp$NS_judge_10_19,.75)),2)




##  percentiles reported at bottom of Figure 4
# 1933-45 share
round(mean(
quantile(cases.temp$NS_judge_1_24,.25),
quantile(cases.temp$NS_judge_2_24,.25),
quantile(cases.temp$NS_judge_3_24,.25),
quantile(cases.temp$NS_judge_4_24,.25),
quantile(cases.temp$NS_judge_5_24,.25),
quantile(cases.temp$NS_judge_6_24,.25),
quantile(cases.temp$NS_judge_7_24,.25),
quantile(cases.temp$NS_judge_8_24,.25),
quantile(cases.temp$NS_judge_9_24,.25),
quantile(cases.temp$NS_judge_10_24,.25)),2)

round(mean(
quantile(cases.temp$NS_judge_1_24,.75),
quantile(cases.temp$NS_judge_2_24,.75),
quantile(cases.temp$NS_judge_3_24,.75),
quantile(cases.temp$NS_judge_4_24,.75),
quantile(cases.temp$NS_judge_5_24,.75),
quantile(cases.temp$NS_judge_6_24,.75),
quantile(cases.temp$NS_judge_7_24,.75),
quantile(cases.temp$NS_judge_8_24,.75),
quantile(cases.temp$NS_judge_9_24,.75),
quantile(cases.temp$NS_judge_10_24,.75)),2)


# pre-1933 share
# 1933-45 share
round(mean(
quantile(cases.temp$NS_judge_1_23,.25),
quantile(cases.temp$NS_judge_2_23,.25),
quantile(cases.temp$NS_judge_3_23,.25),
quantile(cases.temp$NS_judge_4_23,.25),
quantile(cases.temp$NS_judge_5_23,.25),
quantile(cases.temp$NS_judge_6_23,.25),
quantile(cases.temp$NS_judge_7_23,.25),
quantile(cases.temp$NS_judge_8_23,.25),
quantile(cases.temp$NS_judge_9_23,.25),
quantile(cases.temp$NS_judge_10_23,.25)),2)

round(mean(
quantile(cases.temp$NS_judge_1_23,.75),
quantile(cases.temp$NS_judge_2_23,.75),
quantile(cases.temp$NS_judge_3_23,.75),
quantile(cases.temp$NS_judge_4_23,.75),
quantile(cases.temp$NS_judge_5_23,.75),
quantile(cases.temp$NS_judge_6_23,.75),
quantile(cases.temp$NS_judge_7_23,.75),
quantile(cases.temp$NS_judge_8_23,.75),
quantile(cases.temp$NS_judge_9_23,.75),
quantile(cases.temp$NS_judge_10_23,.75)),2)



sessionInfo()





.
